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ABSTRACT 


i/ 

The  rate  of  convergence  of  the  finite  element  method  is  a  function  of 
the  strategy  by  which  the  number  of  degrees  of  freedom  are  increased. 
Alternative  strategies  are  examined  in  the  light  of  recent  theoretical 
results  and  computational  experience. 
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Introduction 


In  finite  element  analysis  convergence  can  be  achieved  in  several 
different  ways,  nevertheless  it  is  useful  to  distinguish  among  three  basic 
modes  of  convergence:  (1)  The  basis  functions  for  each  finite  element  can 
be  fixed  and  the  diameter  of  the  largest  element,  hmny>  allowed  to  approach 
zero.  This  mode  is  called  h- convergence  and  its  computer  implementation 
the  h-version  of  the  finite  element  method.  (2)  The  finite  element  mesh 
can  be  fixed  and  the  minimal  order  of  (polynomial)  basis  functions,  pm1n> 
allowed  to  approach  infinity.  This  mode  is  called  p-convergence  and  its 
computer  implementation  the  p-version  of  the  finite  element  method. 

(3)  Mesh  refinement  can  be  combined  with  increments  in  the  order  of  poly¬ 
nomial  basis  functions.  There  are  many  possible  variations  within  each 
mode  and  there  are  other  convergence  processes  as  well.  For  example,  we 
may  concurrently  decrease  h  x  and  allow  Poisson's  ratio  to  approach  the 
value  of  1/2  to  obtain  the  limiting  case  for  incompressible  solids. 

The  fact  that  convergence  occurs  has  been  the  basis  for  justification 
of  the  finite  element  method  but,  as  far  as  state  of  the  art  finite  element 
analysis  is  concerned,  convergence  is  not  actually  attempted  in  the  compu¬ 
tational  process.  Analysts  are  generally  concerned  with  some  specific 
finite  element  mesh  and  a  corresponding  fixed  set  of  basis  functions.  The 
question  of  whether  the  choice  of  mesh  and  basis  functions  is  adequate  for 
the  purposes  of  an  analysis  is  not  addressed  directly.  Rather,  the  analyst 
relies  on  his  judgement  and  experience  to  ensure  that  the  mesh  is  suffi¬ 
ciently  fine  or  the  polynomial  orders  are  sufficiently  high  so  that  the 
error  of  analysis  is  small.  In  other  words,  finite  element  solutions  are 
intended  to  be  in  the  asymptotic  range  either  with  respect  to  hmav  -*•  0  or 
p^in  -*•  00 •  Because  the  analysts'  judgement  is  not  always  reliable,  there  is 


a  growing  interest  in  adaptivity.  Adaptivity  is  a  procedure  for  efficient 
reduction  of  error  on  the  basis  of  data  already  computed. 

The  most  efficient  error  reduction  technique  is  that  for  which  the 
path  on  the  error  versus  cost  diagram  is  the  steepest.  If  we  simplify  the 
problem  and  assume  that  cost  is  a  simple  ascending  function  of  the  number 
of  degrees  of  freedom  then  the  most  efficient  error  reduction  process  is 
the  convergence  mode  characterized  by  the  highest  rate  of  convergence.  In 
reality  the  cost  depends  on  other  factors  as  well,  some  of  which  are  dif¬ 
ficult  to  quantify,  nevertheless  it  is  useful  to  establish  asymptotic  rates 
of  convergence  for  the  various  modes. 

It  is  possible  to  make  statements  about  asymptotic  rates  of  convergence 
a  priori,  i.e.  without  actually  performing  the  computations,  in  terms  of  a 
property  of  the  approximated  function,  called  "smoothness".  In  this  paper 
we  present  the  definition  of  smoothness  of  functions;  summarize  the  avail¬ 
able  theoretical  knowledge  concerning  the  asymptotic  rates  of  the  various 
modes  of  convergence  and  present  example  problems  from  two-dimensional 
elasticity.  Finally  we  present  some  general  conclusions  concerning  expected 
relative  computational  efficiencies  and  aspects  of  reliability  of  the 
various  modes  of  convergence. 
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2.  Basic  Notation  and  Preliminaries 

Throughout  the  paper,  a  point  in  the  plane  will  be  denoted  by 
x  5  (x1>x2).  n  will  represent  a  polygonal  domain  with  vertices  and 
boundary  T  which  is  the  union  of  (straight)  line  segments  The 

internal  angles  will  be  represented  by  8^.  This  notation  is  shown  in 
Fig.  1. 


We  shall  be  concerned  with  functions  defined  on  Q.  It  will  be 

necessary  to  classify  these  functions  with  respect  to  their  smoothness. 

v 

In  particular,  we  shall  say  that  u  e  H  (ft) ,  k  >  0  integer,  when: 
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where  dx  =  dx^^.  The  Index  k  is  a  measure  of  the  smoothness  of  u  in 
that  it  indicates  how  many  square  integrable  derivatives  u  has.  We 
shall  say  that  a  function  is  smooth  when  k  has  a  high  value. 

As  an  example,  let  us  consider  the  case  of  linear  elastic  fracture 
mechanics,  mode  I:  The  displacement  vector  components  in  plane  strain 
are  [1]: 


ui  "  g1  cos(f)  [1-2'’+sin2  f ]  +  °<r3/2> 

u2  "  G1  ( 3in(f)  [2-2v-cos2  |]  +  °(r3/2> 


2  G 


in  which  r  and  8  are  polar  coordinates,  centered  on  the  crack  tip,  G  and 
v  are  material  constants  and  is  the  stress  intensity  factor.  In  this 
case  both  u^  and  U2  belong  to  H^(ft).  Functions  and  ^  are  members  of 
a  class  of  functions  which  have  great  importance  in  finite  element 
analysis.  The  general  form  of  the  class  is: 


v  -  Re  [rY(logr)6f (0)]  Re[y]  =  a  >  0,  6  >  0 


in  which  f(6)  is  a  smooth  function  of  6.  The  origin  of  the  polar 
coordinates  is  typically  located  at  a  vertex  of  the  domain  ft. 

|l£ 

The  definition  of  the  space  of  functions  H  (ft)  in  eq.  (1)  is  for  k 
integer  only,  nevertheless  it  is  possible  to  generalize  the  notion  of 

lr 

H  (ft)  to  k  >  0  fractional.  See,  for  example,  [2]. 


I 
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s 

-The  function  v,  defined  by  eq.  (3),  belongs  to  all  spaces  H  (ft) , 

s  <  a+1.  In  general,  v  does  not  however  belong  to  H  A(Q).  Thus, 

3  /  2"£ 

u^  and  u^  defined  in  eq.  (2)  belongs  to  H  (ft),  e  >  0  arbitrarily. 

The  foregoing  definitions  provide  for  establishing  two  fundamental 

notions:  If  a  function  belongs  to  H  (ft)  and  the  highest  derivative  in 

the  strain  energy  expression  is  not  greater  than  k,  then  the  function 

has  finite  strain  energy.  Secondly,  the  fractional  index  k  allows  us  to 

quantify  the  idea  of  "smoothness"  of  functions.  For  example  if  ft  contains 

1/3  s 

the  point  r  *  0  then  w^  *  r  cos6  belongs  to  all  spaces  H  (ft),  s  <  4/3, 

1  /  3  s 

and  is  smoother  than  w^  ■  r  cos6  which  belongs  to  all  spaces  H  (ft) , 

s  <  6/5. 

* 

Let  T  =  U  y.,  j  *  l,2,...,m  be  the  union  of  some  sides  of  the 

j  J  k  * 

polygon  ft.  Then  if  u  belongs  to  H  (ft),  k  >  1,  and  u  *  0  on  T  ,  then  we 

write:  u  e  H  ^(ft).  In  the  case  of  vector  functions  u  »  (u^,^), 

u  e  Hk(ft)  means  that  ui  e  Hk(ft),  i  *  1,2. 

The  model  problems  to  be  discussed  have  been  taken  from  two-dimensional 

elasticity,  assuming  plane  strain  conditions.  E  and  v  refer,  respectively, 

to  the  modulus  of  elasticity  and  Poisson's  ratio,  (0<v<l/2) .  Function 

u  ■  (u^,u£)  is  the  displacement  vector  function.  The  strain  energy 

function  is  defined  by: 
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Let  us  now  define  on  r-r  (the  part  of  the  boundary  over  which  no 

displacements  are  prescribed)  the  function  T  =  (t^.t,,),  representing 

* 

prescribed  tractions,  and  on  r  the  function  4  =  (</?,<£.),  representing 

f  2 

prescribed  displacements.  We  assume  that  J  4  t^  ds  <  «  and  there 


«  and  there 


*  1  *  *  1  1  it 

exists  <P  £  H  (ft)  and  ^  ^  on  r  ,  i-1,2.  If  r  j*  t>,  then  we  have  the 

well  known  uniqueness  theorem: 


Theorem  1.  There  exists  exactly  one  function  u  =  (u^,u2>»  u^  e  n(fl) 
such  that: 

* 

(a)  u.  »  (p.  on  T 

l  i 

(b)  u  minimizes  the  potential  energy  functional 


it(u)  -  W(u)  -  /  t  u  d 

i=i  •  * 


among  all  vector  functions  u  *  (u^.u^),  ui  e  h  (ft)  such  that  u^  *  on 
* 

T  .  The  proof  of  this  theorem  is  given  several  textbooks,  see  for 
example,  [3].  If  T  *  (J  and  T  satisfies  the  usual  equilibrium  conditions 
then  u  exists  and  is  determined  by  theorem  1,  except  for  arbitrary  rigid 
body  motion. 

★ 

Assuming  that  T  is  smooth  on  every  side  y  e  r-r  and  $  is  smooth 
* 

on  every  side  y^  e  r  ,  the  solution  u  that  minimizes  tt(u)  in  theorem  1 

is  smooth  on  ft  except  for  possible  singular  behavior  at  the  corners  A^ 

of  ft.  Such  singular  behavior  is  characterized  by  the  coefficients  a  and 

6  in  eq.  (3)  which  depend  on  the  angle  9  and  whether  the  two  sides 

*  * 

adjacent  to  belong  to  r  or  r-r  .  More  detailed  discussion  of 
this  point  is  given  in  [4,5,6].  The  important  observation  is  that  the 


-7- 


at 


space  H  (ft)  to  which  the  solution  u  belongs  is  determined  primarily  by  the 
vertex  angles  of  ft  and  the  boundary  conditions  imposed  on  ft,  provided 
that  the  loading  and  prescribed  displacements  are  smooth. 

The  strain  energy  W(u)  defined  in  eq.  (4)  has  finite  value  for  any 
u  £  H^ft).  Furthermore,  /  W(u)  has  all  of  the  properties  of  a  norm  if 
u  e  H^. (ft)  and  T  0  0.  In  the  following  the  norm  | |u| |  »  /W(u)  plays 

r  E 

essential  role  in  measuring  the  error  of  finite  element  approximation. 

In  fact  | | u | |_  is  equivalent  to  the  norm  | |uj |  ,  l.e.  there  exist  C, , 

H  (ft)  1 

C^,  independent  of  u,  such  that: 


CJMI  1  <  ||u||E  1  C  ||u|| 

H  (ft)  H  (ft) 
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3.  The  Finite  Element  Method 

We  shall  consider  two  different  families  of  meshes:  The  first 

family,  denoted  by  Y,  is  the  family  of  the  usual  triangularizations  with 

bounded  aspect  ratio.  A  specific  triangularization  shall  be  denoted  by 

x.  Thus:  t  e  y.  A  specific  triangle  shall  be  denoted  by  x^.  Thus: 

e  t  e  y  The  largest  side  of  x^  is  denoted  by  lu  and  the  maximum  of 

h . :  h  (t)  «  max  h . ;  the  minimum  of  h. :  h  .  (x)  •  min  h.. 
l  max  ,  i  i  mm  l 


Definition:  The  family  y  is  quasiuniform  if  for  any  x  e  y  we  have: 


h 

max 

h  . 
mm 


(t) 

(xj 


< 


K  <  ® 


The  other  family  to  be  considered  is  a  family  of  square  meshes  G 
(with  possible  refinement).  In  this  case  we  shall  assume  that  all  of  the 
sides  of  ft  are  parallel  with  the  coordinate  axes.  Although  this  assumption 
may  seem  to  be  overly  restrictive  from  the  practical  point  of  view,  in 
reality  it  is  not:  Any  plane  domain  can  be  subdivided  into  rectangular 
figures  with  curved  (or  straight)  sides  such  that  each  figure  can  be 
mapped  onto  the  unit  square  by  smooth  variable  transformations. 

An  example  of  the  square  meshes  is  given  in  fig.  2.  The  nodal  points 
of  the  mesh  which  are  either  lying  on  T  or  are  common  vertices  of  four 
squares  are  called  regular,  the  others  are  called  irregular.  The  term: 
"quasiuniform  mesh"  is  applied  to  the  square  meshes  in  the  same  way 
as  before,  i.e.  it  means  that  the  ratio  of  the  largest  element  diagonal 
to  the  smallest  is  bounded  as  the  size  of  the  largest  element  approaches 


zero. 


Fig.  2 

Square  mesh  with  refinement 

Over  each  finite  element  we  approximate  u  =  (u^.u^),  u  e  H^(fi),  with 

polynomials  of  order  p  and  we  assume  that  u  e  H1^).  The  polynomials 

* 

assume  the  value  of  $  *  °n  r  *  The  P0^^0111^^  order  p  may  vary 

from  element  to  element.  In  the  following  we  shall  assume  also  that  the 
interelement  continuity  requirements  are  enforced  exactly  and  minimally, 
i.e.  overconformity  and  nonconformity  are  avoided.  (The  assumption  that 
u  e  H^(n)  guarantees  conformity).  The  number  of  degrees  of  freedom, 
after  enforcement  of  the  principal  boundary  conditions,  is  denoted  by  N. 

We  shall  denote  the  set  of  functions  u  e  H^(ft)  which  are  polynomials 
of  degree  at  most  p^  on  every  e  t  and  satisfy  the  principal  boundary 

•k 

conditions  ♦  on  F  by  =  M^(t ,p,$) .  (The  polynomial  order  may  vary 
from  element  to  element).  Similarlv,  we  shall  define  the  set  of  functions 
u,  which  are  bilinear  (i.e.  products  of  linear  polynomials  in  x^  and 

"ft 

x2)  on  every  element  of  a  square  mesh  and  u  »  ♦  on  F  by  Np^. 


The  finite  element  solution  e  Mj,^  (respectively  e  N^)  is 
the  function  that  minimizes  the  energy  functional  given  by  eq.  (5)  over 
the  set  (respectively  . 

Let  u  be  the  exact  solution  by  theorem  1.  Then  it  is  not  difficult 
to  prove  that: 

llU-UFE^2*  |^(u)-«(uFE)  j  -  |W(u)-W(upE)  |  (6) 

E 

i.e.  the  energy  of  the  error  of  the  finite  element  solution  equals  the 
error  of  the  energy  [ 7 ] . 
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4.  The  h  and  p  Versions  of  the  Finite  Element  Method 

In  the  classical  finite  element  method  the  polynomial  degree  is 

fixed,  usually  at  some  low  p  value  (p=l,2,  or  3)  and  the  solution  is 

constructed  in  M_.(t,p,$)  such  that  h  is  small.  Increased  accuracy 
l  *  max 

is  achieved  by  mesh  refinement  (h  -*■  0) .  This  strategy  for  increasing 

max 

the  number  of  degrees  of  freedom  is  called  the  h-version  of  the  finite 
element  method. 

We  can  also  fix  the  mesh  t  (with  relatively  small  number  of 
triangles)  and  increase  p,  either  uniformly  or  selectively.  This 
strategy  for  increasing  the  number  of  degrees  of  freedom  is  called  the 
p-version  of  the  finite  element  method. 

The  h  and  p  versions  can  be  viewed  as  special  cases  of  the  finite 
element  method  which  allows  increasing  the  number  of  degrees  of  freedom 
by  concurrent  refinement  of  the  mesh  and  increases  in  p. 
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5.  Model  Problems 

We  now  define  the  model  problems  used  for  illustrating  the  convergence 
properties  of  the  various  strategies  used  for  increasing  the  number  of 
degrees  of  freedom  in  the  finite  element  method. 

Problem  1;  Square  domain  under  plane  strain  conditions  subjected  to 
imposed  shear  displacement,  ft  is  the  square,  as  shown  in  Fig.  3. 


Fig.  3 

Square  domain,  example  problem  1. 


r  “  Y1  U  Y3 


!+l  on  y, 

-1  on  y. 


*  0  on  U  Y  3 
T  -  (0,0)  on  y 2  (J 


We  shall  be  concerned  with 
the  case  of  E  »  1,  v  ■  0.3 
and  v  *  0.4999.  Because 
and  x^  are  axes  of  anti¬ 
symmetry,  we  shall  compute 
the  finite  element  solution 


Upg  for  the  right  upper  quarter  of  the  domain  only.  Consequently  the 
number  of  degrees  of  freedom  shall  be  related  to  this  quarter  only. 
The  estimated  exact  strain  energies  for  the  quarter  domain  are: 


W  ■  0.130680  for  v  -  0.3  and  W  ■  0.127035  for  v  ■  0.4999. 
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The  solution  u  has  singular  behavior  in  the  neighborhood  of  all 
vertices  of  the  form  given  in  eq.  (3)  (with  the  origin  of  the  polar 
coordinate  system  at  vertex  A^) .  We  observe  that  u  e  H  (Q) ,  where  k 
depends  on  Poisson's  ratio,  as  shown  in  fig.  4.  The  data  in  fig.  4  were 
computed  on  the  basis  of  reference  4. 


Fig.  4 

The  smoothness  parameter  k  =  1+a  as  a  function 
of  Poisson's  ratio  in  example  problem  1 

Problem  2.  The  edge  cracked  square  panel  under  plane  strain  conditions 
subjected  to  uniform  tension  as  shown  in  Fig.  5.  In  this  case  r  ■  t. 
The  boundary  conditions  are: 

T  “  (0,0)  on  T  -  y1  |j 

T  -  (0,1)  on 

T  -  (0,-1)  on  y6 
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We  shall  be  concerned  the  case  of  E  *  1,  v  ■  0.3  and  v  ■  0.4999.  Because 
x^  and  are  axes  of  symmetry,  we  shall  compute  the  finite  element 
solution  Upg  for  the  right  upper  quarter  of  the  domain  only  and  the 
number  of  degrees  of  freedom  shall  be  related  to  this  domain  only.  The 
estimated  exact  strain  energies  for  the  quarter  domain  are:  W  -  0.73422 
for  v  *  0.3  and  W  *  0.60525  for  v  *  0.4999. 


Fig.  5 

Edge  cracked  square  panel 

The  solution  u  has  singular  behavior  at  the  crack  tip  (vertices  A^,  and 

■i/2-c 

Ag)  of  the  form  given  by  eq's.  (2a, b).  Thus  u  e  H  '  (0).  Unlike  in 

example  problem  1,  the  smoothness  parameter  k  is  independent  of  Poisson's 


ratio. 
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6.  Asymptotic  Error  Analysis 

When  comparing  alternative  approaches  in  the  finite  element  method, 
the  error  vs.  cost  relationship  is  of  interest,  with  the  error  measured 
in  terms  of  displacements,  stresses,  stress  resultants  at  specific 
points,  or  in  terms  of  energy.  The  error-cost  relationship  is  often 
simplified  by  making  the  assumption  that  the  cost  is  some  simple  function 
of  the  number  of  degrees  of  freedom.  We  shall  follow  the  same  course 
and  use  the  error  vs.  number  of  degrees  of  freedom  as  our  basis  for 
comparison.  Thus  a  given  solution  will  be  represented  by  a  point  on  the 
error  vs.  N  diagram.  We  shall  be  concerned  with  the  rate  of  change  of 
the  error  in  energy  with  respect  to  N.  N  can  be  increased  in  various 
ways:  uniform  or  quasiuniform  mesh  refinement;  non-quasiuniform  mesh 
refinement,  uniform  or  non-uniform  change  in  polynomial  order  on  a  fixed 
mesh,  etc.  These  can  be  viewed  as  various  "extensions"  of  the  original 
solution.  We  shall  say  that  an  extension  is  asymptotic  or  pre-asymptotic 
depending  on  whether  the  error  vs.  N  relationship  is  governed  by  one  or 
more  parameters.  When  N  is  sufficiently  large  then  the  error  can  be 
characterized  well  by  a  function  t|>(N)  •  Function  ip  depends  in  general 
on  the  smoothness  of  the  approximated  function.  Thus  the  error  can 
be  written  as: 

|  | e j  |£  -  *<H)+0[*<N)] 

ip  (N)  =  CN~a 


with: 


where  a  (a>0)  is  Che  governing  parameter  and  o[i|>(N)]  represents  terms 
that  approach  zero  as  N  +  •  faster  than  <HN).  In  order  to  ensure  that 
the  various  extensions  are  compared  on  the  same  basis,  we  shall  use  a  as 
our  basis  for  comparison.  Thus  we  shall  make  the  assumption  that  each 
extension  is  in  the  asymptotic  range.  This  is  important  from  the  point 
of  view  of  reliability:  In  the  asymptotic  range  the  error  is  not  sig¬ 
nificantly  influenced  by  optional  input  parameters  (such  as  the  kind  of 
mesh  divisions  used)  and  the  error  is  generally  small. 


I 


7.  The  h-version 


We  shall  now  summarize  theorems  concerning  the  asymptotic  rate  of 
convergence  In  the  h-verslon  of  the  finite  element  method. 

Theorem  2.  Let  the  exact  solution  u  belong  in  the  space  H  (£2) .  Then  for 
a  family  of  quasiuniform  meshes  y  and  the  space  of  polynomial  approxi¬ 
mating  functions  Mr  .  (t  ,  p,4>)  >  in  which  p  is  fixed  and  uniform,  we  have: 
r* 

Mu-UpJI  <  C(n,k,p,Y)S“  1/2  mln(k‘1»p)|  |u|  I  (7) 

e  "  h  (n) 

In  eq.  (7)  the  constant  C  does  not  depend  on  u  or  N  [2,8].  The  absolute 
value  of  the  exponent  of  N  is  called  the  asymptotic  rate  of  convergence 
or  simply  rate  of  convergence. 

An  inverse  theorem  exists  also,  which  can  be  summarized  as  follows: 

If  we  were  able  to  observe  the  asymptotic  rate  of  convergence  for  a  given 
problem  for  any  quasiuniform  triangular  mesh  and  fixed  p  (assuming  complete 
polynomials)  to  be  1/2  a  (i.e.  a  *  y  a) : 

1 !u-ufeI I  i  ctf1/2a 

ra  E 

then  we  would  be  able  to  make  the  following  statements  about  the  exact 
solution  u: 

1)  If  1  <  a  <  p  then  u  e  H1+a""£  (fi) ,  e  >  0  arbitrarily. 

2)  If  a  >  p  then  u  is  a  polynomial. 

3)  If  a  »  p  then  u  e 

For  analysis  of  the  inverse  theorem  we  refer  to  [9,10]. 

Theorem  2  and  its  inverse  (summarized  in  the  preceding  statements) 
show  that  the  asymptotic  rate  of  convergence  is  completely  characterized 
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by  the  smoothness  of  u.  In  most  problems  of  practical  importance  k  is 
between  1.5  and  2,  and  u  is  not  a  polynomial,  hence  the  asymptotic  rate 
of  convergence  is  governed  by  k.  Theorem  2  indicates  that  changing  p 
when  p  >  k-1  will  not  change  the  asymptotic  rate  of  convergence,  but 
will  affect  the  accuracy  of  the  solution  because  C  depends  on  p. 

We  shall  now  illustrate  theorem  2  on  problem  1  with  v  ■  0.3.  The 
uniform  mesh  for  the  quarter  domain  (with  h  variable)  is  shown  in  Fig.  6. 


Fig.  6 

Uniform  triangulation  (Example  problems  1  and  2) . 


The  computations  for  the  uniform  meshes  were  performed  by  means  of 
COMET-X,  a  computer  program  developed  at  Washington  University.  COMET-X 
implements  the  p-version  of  the  finite  element  method  through  the  use  of 
exactly  and  minimally  conforming  hierarchic  finite  elements.  The  range 
of  polynomials  permitted  by  COMET-X  is  from  1  to  8  [11,12], 
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Adaptively  constructed  non-quasiuniform  meshes  are  shown  in  Fig.  7. 
These  meshes  were  generated  and  the  computations  performed  by  FEARS 
(Finite  Element  Adaptive  Research  Solver),  a  computer  program  developed 
at  the  University  of  Maryland.  FEARS  gives  a  reliable  estimate  of 
error  of  the  finite  element  solution  [13]. 

The  relative  error  is  plotted  against  N  on  a  log-log  scale  in  Fig.  8 
for  various  p  values.  The  slope  of  these  curves  for  large  N  is  the 
asymptotic  rate  of  convergence.  We  note  that  in  Fig.  8  (and  in  the 
other  figures  illustrating  the  error  vs.  N  relationship)  we  measure  the 
error  by  the  square  of  the  norm  ||*||g.  See  eq.  (6). 

It  is  seen  that  for  uniform  mesh  refinement  the  error  depends  on 
p  but  the  slope  is  independent  of  p.  With  k  =  1+a  =  1.76  from  Fig.  4, 
the  slope  of  0.76  is  as  predicted  by  theorem  2. 

It  has  been  proven  that  if  the  approximated  function  u  is  of  the 
functional  form  of  eq.  (3)  then  there  exists  a  sequence  of  non-quasiunlform 
meshes  such  that  the  rate  of  convergence  depends  only  on  p,  not  on  y  and  <5 
as  would  be  the  case  if  uniform  or  quasiuniform  meshes  were  used.  The 
proof  is  available  in  [8,10].  The  sequence  of  meshes  having  this  character 
can  be  constructed  a  priori  or  adaptively.  The  meshes  constructed  a 
priori  cause  higher  approximation  error  than  adaptively  constructed 
ones ,  however . 

As  seen  in  Fig.  8,  the  error  vs.  N  curve  for  adaptively  constructed 
meshes  approaches  the  slope  of  -1  for  log-log  scale,  which  is  the  maximum 
rate  of  convergence  possible  for  p  ■  1. 


m  meshes 
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Fig.  9 


Example  problem  1:  Relative  error  in  energy  vs.  number 
of  degrees  of  freedom.  H-version,  Poisson's  ratio:  0.4999. 
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Eq.  7  is  valid  for  all  Poisson's  ratios  0  <  v  <  1/2,  with  the 
constant  C  possibly  dependent  on  v.  Fig.  9  shows  the  behavior  for 
v  ■  0.4999.  We  see  that  the  rate  of  convergence  approaches  the  value 
of  k-1  (in  this  case  0.69)  in  accordance  with  theorem  2  but  for  p  •  1 
the  asymptotic  range  is  entered  at  high  N  values  (in  this  case  beyond 
the  range  of  plotted  values  and  probably  beyond  the  round-off  limitations 
of  digital  computers) . 

The  nearly  degenerate  case  of  v  a  1/2  deserves  special  consideration 
because  two  parameters  are  involved;  1/2  -  v  and  h.  The  asymptotic 
theory  is  applicable  only  when  h  is  small  with  respect  to  1/2  -  v.  It 
is  well  known  that  elements  with  p  =  1  perform  poorly  when  v  a  1/2  and 
various  special  approaches,  such  as  reduced  integration,  have  been 
proposed.  The  results  shown  in  Fig.  9  indicate  that  mesh  refinement 
will  not  reduce  the  error  when  p  =  1,  on  the  other  hand  use  of  p  >  3 
essentially  eliminates  this  difficulty.  Theoretical  analysis  of  this 
effect  is  not  yet  available. 


8.  The  p-version 

We  shall  now  review  the  basic  properties  of  the  p-version. 


Definition:  Let  p  denote  the  largest  polynomial  order  of  the  basis 
-  max 

functions  over  all  finite  elements  and  let  p  ,  denote  the  smallest.  A 

mm 

sequence  of  p-distributions  is  quasiuniform  if: 


-  <  K  <  <*> 

pmin 

Theorem  3.  Let  the  exact  solution  u  belong  in  the  space  H  (ft) .  Then 
for  the  space  of  exactly  and  minimally  conforming  polynomial  approximating 
functions  Mj,A(t,p,$)  in  which  the  triangular  mesh  T  is  fixed  and  the  sequence 
of  p-distributions  is  quasiuniform,  for  any  e  >  0  we  have: 

Mu-Up- 1  |  <  C(ft,k,T,K,e)N'1/2(k"1)+£M“M  k  (8) 

Theorem  3  is  similar  to  theorem  2,  however  here  C  is  independent  of 

p.  The  proof  is  given  in  [14] .  Importantly,  the  inverse  theorem,  also 

—1/  2rt 

given  in  [14],  states  that  if  | |u-UpE| |  <  CN  under  the  conditions  of 

1+a-e  *  * 

theorem  3,  then:  (a)  u  e  H  (ft  )  where  ft  is  a  subdomain  of  ft  not 

1  l  /  9 

containing  any  of  the  finite  element  boundaries;  (b)  u  e  H  E(ft). 

This  is  significantly  different  from  the  inverse  of  theorem  2.  In 

particular,  if  the  singular  behavior  of  u  is  confined  to  the  boundaries 

* 

of  t  (i.e.  the  singularity  is  not  in  ft  )  then  the  rate  of  convergence  of 
the  p  version  is  twice  the  rate  of  convergence  of  the  h-version,  provided 
that  p  >  k-1  in  the  h-version.  In  other  words,  the  p-version  can  "absorb" 
singular  behavior  at  the  boundaries  of  finite  elements. 


A  very  important  case  is  when  the  singularity  is  at  the  corner  of 
the  domain  (and  therefore  at  the  vertex  of  one  or  more  elements).  Then 
we  can  prove  the  following  theorem: 


Theorem  4.  Let  u  *  u^-H^,  u^  having  the  functional  form  given  in  eq.  (3) 
with  the  origin  at  the  vertex  of  the  domain  and  u^  e  H  (ft).  Then  for 
triangular  meshes: 


<  C(T,k,ate){"N~a+e|  |u  |  j 
E  “  L  1  H°(ft) 


+  „-l/2(k-l)+e 


Hk(fl) 


] 


(9) 


A  conjectural  statement  concerning  the  existence  of  theorem  4  and  a 
numerical  demonstration  of  this  theorem  was  given  in  [15].  The  proof  is 
given  in  [14]. 

In  summary,  the  p-version  cannot  have  lower  rate  of  convergence 
than  the  h  version  based  on  quasiuniform  meshes.  The  rate  of  convergence 
of  the  p-version  is  twice  that  of  the  h-version  when  the  singularity  is 
at  element  boundaries  and  quasiuniform  meshes  are  used.  On  the  other 
hand  the  h-version,  with  the  use  of  optimally  refined  meshes  (which  are 
not  quasiuniform)  and  sufficiently  high  p  can  have  higher  rate  of 
convergence  than  the  p-version. 

Fig.  10  illustrates  that  the  rate  of  convergence  of  the  p-version 
is  twice  the  rate  of  convergence  of  the  h-version  based  on  uniform  mesh 
refinement  in  the  case  of  example  problem  2.  The  fact  that  the  asymptotic 
range  is  entered  at  low  p  values  is  noted. 

In  Fig.  10  we  also  show  the  performance  of  adaptively  constructed  square 


meshes  for  p  *  1.  Once  again  we  see  that  the  error  vs.  N  curve  approaches 
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Fig.  10 

Example  problem  2:  Relative  error  in  strain  energy  vs.  number 
of  degrees  of  freedom.  Comparison  of  the  h  and  p  versions. 
Poisson's  ratio:  0.3 
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Adaptively  constructed  non-quasiuniform  meshes 
for  example  problem  2. 
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Fig.  12 


Example  problem  2:  Relative  error  in  strain  energy  vs.  number  of 
*  degrees  of  freedom.  Comparison  of  the  h-  and  p-versions. 

Poisson’s  ratio:  0.4999. 


RELATIVE  ERROR  (PERCENT) 
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che  slope  of  -1  on  log-log  scale  for  large  N,  which  is  che  largest 
asymptotic  rate  possible  for  p  ■  1.  The  larger  pre-asymptotic  rate  is 
due  to  the  fact  that  at  the  beginning  a  very  small  number  of  additional 
elements  at  the  crack  tip  increase  the  accuracy  significantly.  Subse¬ 
quently  the  domain  must  be  refined  away  from  the  singularity,  as  shown 
in  Fig.  11,  which  causes  the  rate  of  convergence  to  become  slower. 

Fig.  12  illustrates  that  the  point  of  entry  into  the  asymptotic 
range  is  not  affected  in  the  p-version  to  a  significant  degree  when 
Poisson's  ratio  is  close  to  0.5,  but  is  strongly  affected  in  the 
h-version  when  p  *  1.  As  in  example  1,  mesh  refinement  is  not  effective 
for  error  reduction. 
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9.  Rate  of  Convergence  when  pmav  is  Increased  with  Concurrent 

Non-quasiuniform  Mesh  Refinement. 

In  this  section  we  analyze  some  convergence  properties  of  the 
finite  element  method  under  the  condition  that  mesh  refinement  is 
accompanied  by  increases  in  p.  First  we  quote  the  following  basic 
theorem: 

Ir 

Theorem  3:  Let  u  e  H  (Q) ;  let  y  be  a  family  of  quasiuniform  triangular 
meshes  and  let  u^  be  the  finite  element  solution  based  on  the  space  of 
piecewise  polynomial  functions  (t,p,4>)  with  t  e  y;  P  quasiuniformly 
distributed,  and  Pmin  >  k-1.  Then,  for  any  e  >  0 

-ki+e 

||u-u_  I!  <  C(fl,k,e)N  2  | I U I  I  (10) 

E  "  IT(ft) 

The  proof  of  this  theorem  is  given  in  [16].  The  significance  of  this 
theorem  is  that  it  joins  theorems  3  and  4,  making  the  constant  C  inde¬ 
pendent  of  both  p  and  t.  An  inverse  theorem  exists  but  Is  not  known 
as  yet. 

In  Sections  7  and  8  we  demonstrated  that  it  is  possible  to  generate 
proper  sequences  of  mesh  such  that  the  rate  of  convergence  is  independent 
of  the  singularity  but  depends  on  the  p-distribution.  In  the  case  of 
corner  singularities  (of  the  form  given  by  eq.  3)  it  can  be  shown  (see 
[16])  that  properly  refined  sequences  of  mesh  combined  with  suitable 
sequences  of  p-distribution  result  in  the  error  bound: 

I  1 u-uFE l |  1  C(B)N'e 

E 


(11) 
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with  0  arbitrarily  large.  It  is  also  known  that  under  certain  conditions 
the  bound  is  exponential  [16]. 

We  shall  demonstrate  on  example  problem  1  that  6  can  be  arbitrarily 
large.  The  sequence  of  meshes  and  corresponding  p-dlstrlbutions  are 
shown  in  fig.  13.  Each  refinement  reduces  the  size  of  the  triangles  at 
the  singularity  by  the  factor  (1-p).  Thus  the  corresponding  sides  of 
the  triangles  at  the  singularity  are  in  geometric  progression  with 
common  ratio  (1-p).  We  note  that  this  refinement  is  not  quasiuniform. 

The  p-distribution  is  also  not  quasiuniform:  p  -  1  for  the  two  elements 
at  the  singularity;  p  *  2  for  the  next  group  of  four  elements  then  p 
progressively  increases  by  increments  of  one  for  each  additional  group 
of  four  elements  away  from  the  singularity. 

The  relative  error  in  strain  energy  vs.  N  is  plotted  on  log-log 
scale  in  Figures  14  and  15  for  two  refinements  characterized  by  p  *  0.62 
(the  "golden  rule"  refinement  which  is  shown  in  Fig.  13)  and  p  -  0.90,  a 
much  stronger  refinement.  It  is  seen  that  the  slope  of  the  relative 
error  vs.  N  curve  progressively  increases  with  N,  indicating  that  3  in 
eq.  (11)  is  an  ascending  function  of  N  and  can  be  arbitrarily  large. 

The  results  obtained  with  the  p-version  for  two  elements  are  also 
shown  in  Figures  14  and  15.  The  indications  are  that  within  the  range 
of  accuracy  of  practical  interest  the  p-verslon  is  as  effective  in 
reducing  the  error  as  the  strategy  just  outlined. 

Optimal  combination  of  non-quasiunlform  mesh  refinement  and 
p-distribution  is  a  delicate  matter.  For  very  high  accuracy  one  can 
expect  the  polynomial  degrees  for  elements  at  the  singularity  to  be 
smaller  than  for  the  larger  triangles  away  from  the  singularity.  Never¬ 
theless,  in  the  pre-asymptotic  range  the  optimal  p-distribution  can  be 
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Fig.  13 

Example  problem  1:  Non-quasiuniform  mesh  refinement  and  p-distribution. 
(The  p-values  are  shown  for  each  element) .  The  18-element  mesh  with 
o  -  5  is  not  shown. 
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Example  problem  1:  Relative  error  in  strain  energy  vi 
number  of  degrees  of  freedom  for  non-quasiuniform  mei 
refinement  and  p-distribution.  Poisson's  ratio:  0.3 
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Fig.  15 

Example  problem  1:  Relative  error  in  strain  energy  vs.  number 
of  degrees  of  freedom  for  non-quasiuniform  mesh  refinement  and 
p-distribution.  Poisson's  ratio:  0.4999 


RELATIVE  ERROR  (percent) 
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quite  different.  For  general  theoretical  results  and  some  numerical 
results  in  one  dimension  the  reader  is  referred  to  [16]. 
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10.  Adaptivity 

We  have  seen  that  if  the  number  of  degrees  of  freedom  is  increased 
through  proper  combinations  of  mesh  refinement  and  p-distribution,  the 
error  will  decrease  very  rapidly  with  increasing  N.  The  choice  of 
proper  mesh  refinement  and  p-distribution  depends  on  the  smoothness  of 
the  solution,  however,  and  cannot  in  general  be  determined  a  priori. 

It  is  possible  to  compute  local  measures  of  error,  which  Indicate 
the  contribution  of  each  element  to  the  total  error  of  approximation. 

The  local  error  measures  provide  a  basis  for  establishing  proper  distri¬ 
butions  of  the  degrees  of  freedom. 

A  more  detailed  discussion  of  adaptivity  will  be  presented  in  a 
future  paper. 


11.  Conclusions 


We  can  summarize  Che  main  results  concerning  asymptotic  rates  of 

convergence  in  the  finite  element  method  as  follows: 

1.  The  asymptotic  rate  of  convergence  depends  on  the  smoothness  of  the 
function  to  be  approximated  and  the  order  of  the  polynomial  basis 
functions  (p) .  Smoothness  is  measured  by  the  number  of  square 
integrable  derivatives  (k)  over  the  domain  of  interest,  with  k 
generalized  to  fractional  values. 

2.  In  the  h-version  of  the  finite  element  method  the  rate  of  convergence 
is  the  smaller  of  p  and  k-1  if  uniform  or  quasiuniform  mesh  refine¬ 
ment  is  used.  When  the  singular  behavior  is  caused  by  corners, 
there  is  a  sequence  of  not  quasiuniform  meshes  (called  proper  mesh 
refinement)  for  which  the  rate  of  convergence  is  dependent  only  on  p. 

3.  In  the  p  version  of  the  finite  element  method  the  rate  of  conver¬ 
gence  cannot  be  slower  than  in  the  h-version,  provided  that  quasi¬ 
uniform  mesh  refinement  is  used  in  the  h-version.  When  the  singu¬ 
larity  is  at  the  boundaries  of  finite  elements,  the  rate  of  con¬ 
vergence  of  the  p-version  is  twice  that  of  the  h-version  provided 
that  quasiuniform  mesh  refinement  is  used  and  p  >  k-1  in  the 
h-version. 

4.  It  is  possible  to  design  optimal  sequences  of  meshes  and  p-distribu- 
tlons  for  which  the  rate  of  convergence  in  the  presence  of  singu¬ 
larities  is  arbitrarily  large;  in  fact  the  convergence  can  be 
exponential.  Such  sequences  are  not  quasiuniform  and  depend  on 

the  function  to  be  approximated.  For  this  reason,  the  sequences 
can  be  determined  in  practice  only  by  an  adaptive  approach. 
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For  obvious  practical  reasons,  finite  element  analyses  should  be 
both  efficient  and  reliable.  This  means  that  the  point  of  entry 
into  the  asymptotic  range  should  occur  at  small  values  of  N  and 
should  not  be  sensitive  to  the  input  parameters.  The  p-version 
meets  this  requirement  in  general  better  than  the  h-version.  This 
was  demonstrated  through  the  example  of  nearly  incompressible 
solids,  in  which  the  point  of  entry  into  the  asymptotic  range  was 
not  affected  in  the  p-version  to  an  important  degree  but  was  signi¬ 
ficantly  shifted  in  the  h-version. 

Our  discussion  and  comparisons  were  based  on  error  vs.  number  of 
degrees  of  freedom  relationships.  This  is  a  simplified  treatment 
of  the  more  important  error  vs.  cost  relationship.  Detailed  analysis 
of  marginal  cost  vs.  error  reduction  is  beyond  the  scope  of  this 
paper,  but  we  note  that  reasonably  designed  fixed  meshes,  combined 
with  uniform  or  selective  increases  in  p  provide  the  most  promising 
approach  to  efficient  quality  control  in  finite  element  analysis. 
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